Numerical Evaluation of the “Dual-Kernel 
Counter-flow” Matric Convolution Integral 
that Arises in Discrete/Continuous (D/C) 

Control Theory 

Douglas D. Nixon 
NASA, Marshall Space Flight Center 
Huntsville, Alabama 35812 USA 


,4 Av/rac/- Discrete/Continuous (D/C) control theory is a 
new generalized theory of discrete-time control that expands 
the concept of conventional (exact) discrete-time control to 
create a framework for design and implementation of discrete- 
time control systems that include a continuous-time command 
function generator so that actuator commands need not be 
constant between control decisions, but can be more generally 
defined and implemented as functions that vary with time 
across sample period. Because the plant/control system 
construct contains two linear subsystems arranged in tandem, 
a novel “dual-kernel counter-flow” convolution integral 
appears in the formulation. As part of the D/C system design 
and implementation process, numerical evaluation of that 
integral over the sample period is required. Three 
fundamentally different evaluation methods and associated 
algorithms are derived for the constant-coefficient case. 
Numerical results are matched against three available 
examples that have closed-form solutions. 

I. Introduction 

As introduced in [1], Discrete/Continuous (D/C) control 
theory defines a new generalized approach to design of 
discrete-time control systems that, unlike traditional 
discrete-time control, does not restrict actuator commands to 
be constant over the control decision interval , but allows 
them to vary continuously with time, in a “smart” manner, 
across each sample period between decisions. Because the 
system construct contains two linear subsystems that 
operate in tandem, a novel dual-kernel counter-flow matric 
convolution integral naturally arises in the formulation that, 
in its most general form, is a matrix defined by 

h 

D (t 2 A) = { l 0 p » T ) B p ( T ) C e ( T ) 0 c (A » (!) 

h 

where (A,/, ) are times associated with extremes of the 
sample period, (d> p ,d> c ) are state transition matrices 
associated with the plant and control system, B p (t) is the 
plant input distribution matrix, C c (t) is the control system 
output selection matrix, and the definition and form have 
been taken from [1], 


Specializing (1) to the constant coefficient case yields 

D(T) = l eAp ™BpC c eAclr) dT , (2) 

o 

where the constants A and B are associated with the 
standard state-variables representation of linear system 
dynamics as indicated by 


x = Ax+Bu 

(3) 

6 

H 
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the subscripts p and c, as before, refer to “plant” and 
“control”, and T is the time between control decisions 
(sample period). A block diagram representation of the 
system theoretical construct is shown in Fig. 1 . 

Design and computation of the control gain-matrix 
K shown in Fig. 1 requires evaluation of the convolution 
integral D [1, 2]. The rest of this paper is devoted to 
development of three specific approaches to numerical 
evaluation of D . The first method relies on series 
expansion and truncation for the individual plant and 
control system state transition matrices, and also makes use 
of the state transition matrices to partition the sample 
period so that the length of time for which the series 
expansions must be valid can be less than the sample 
period. A second method, which can be described as a 



Fig. 1. Diagram of Plant/Control System 












mapping-matrix approach, requires simulation of the system 
to generate plant time responses across the sample period as 
a function of control system initial conditions. The 
simulation is run as many times as there are control system 
states. The third method relies on a single series expansion 
and truncation for the total system state transition matrix 
followed by extraction of the convolution as a partitioned 
block of the matrix result. While this method involves 
indirect computation of D from the solution of the total 
system matric differential equation, Johnson [3] describes a 
potential fourth approach that involves identification of a 
matric differential equation whose solution directly 
produces D . 

II. Development of Methods and Algorithms 


A. Method 1: Sample-Period-Partitioned Double Series 
Expansion 

The first method and algorithm are derived by working 
directly with the definition of D given in (2). Both 
exponential functions are replaced with their infinite series 
counterparts so that 


D = e 1 


]\f(~ A P*r 

1 

R C 

00 A ^ m 

V A c r J 

J ^ n\ 

o L n =° 

J 

s 

^1 


dr. (5) 


Exchanging integration with summation, a single element of 
the resulting (double) summation is expressed by, 
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where d mn is any element of D . Now the integration can 
be performed analytically, and the corresponding general 
element of D is given by 
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Finally, the integral D can be expressed as 
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where the integers N and M represent the arbitrarily large 
but finite number of terms included in the series plant and 
(P c of the control system, expansions for the state transition 
matrices ® P of the 

Equation (8) potentially represents an algebraic 
definition of requirements for the detailed design of a 
computer algorithm, which can be implemented in a 


language such as Fortran, C, or Matlab®, to numerically 
evaluate D . However, for reasons of practical experience, 
the algorithm derivation will be taken through another 
major step. The time between control decisions T will be 
partitioned into an arbitrary number of sub-intervals q, and 
the convolution integral D(T) will be expressed as a 
function of D(T / q) and the state transition matrices 
® P (T / q) and ® c (T/q). Because the final symbolic 
expression for D(T) will explicitly include both the 
number of subintervals q and the number of terms N and 
M used in the exponential series expansions, the algorithm 
will become more tunable for specific circumstances and 
much more capable of handling long sample periods. 
Choosing values for the tuning parameters N,M, and q is 

loosely analogous to choosing order in time and length of 
time-step associated with numerical integration of 
differential equations. This form of the algorithm will be 
referred to as sample-period-partitioned (SPP). 

To facilitate derivation of the SPP algorithm, a standard 
compact discrete-time control notation is employed and 
explained below. For constant input u over the sample 
period T , often referred to as sample-and-hold or zero- 
order-hold (ZOH), the exact discrete-time equivalent to 
the continuous-time system represented by (3) is 

x{k + 1) = A{T)x{k) + B(T)u(k) (9) 

for k = 0, 1, 2, ••• , where A is the state-transition matrix 
A(T) = e AT , (10) 

and B is the discrete input distribution matrix 

T 

B{T) = ^e A{T r) Bdr . ( n ) 

o 

The tacit assumption of constant sample period, T(k) = T , 
in (9-11) is a notational convenience that affords no 
significant loss of generality. 

As developed in [1], the D/C generalization of (9) for 
the plant is 

(k + 1) = A p ( T)x p (k) + D(T)x c ( k ) , (12) 

where the value of x c ( k ) is updated from the D/C control 
computations at each sample-time t k = kT based on 
x p ( k ) = x p ( kT) and the D/C control gains associated with 

a particular design. The continuous variation of the D/C 
control system state as it is propagated across the sample 
period is given by 

x c {t) = e Ac(t - m x c {kT)\ kT <t<(k + \)T. (13) 



In the more compact notation of (12), an evaluation of _ JL^ „ — „ _ n 

the control state at the end of the sample period, prior to any D(T) = y A p (T / q)]D(T / q)A 9 ( Ttq ), (23) 

potential alteration (update) by the D/C control logic, is n=1 

given by where, from (8) 


x c (k + X) = A c (T)x c {k). (14) 

Equations (12) and (14) form a basis from which the 
sample period can be partitioned. To this end, an integer q 
is defined such that it represents the number of sample 
period sub-intervals resulting from the partitioning. Then 
(12) applied to the first sub-interval becomes 

x p (k+l/q) = A p (T/q)x p (k) + D(T/q)x c (k) . (15) 

Similarly, (14) becomes 

x c (k + l/q) = A c (T/q)x c (k). (16) 

Continuing forward in time, 


D(T!q) = 

C MTK*YY (-1 TK B P C c A : ( T/q) m+n+l (24) 

m\n\(m + « + 1) 

and it should be pointed out that the notation (T/q) on the 
right-hand side of (24) simply indicates an algebraic 
parameter grouping, and is different from but consistent 
with the functional notation (T/q) used previously and on 
the left-hand side of (24) . The complete sample-period 
partitioned algorithm, defined by (23) and (24), can now be 
applied over a very long control system decision-making 
interval (large T) because the computational process is 
explicitly tunable in terms of both the order M and N of 
the series expansions and the associated time interval T / q . 


x p (k + 2/q) = A p (T/q)x p (k + l/q) 

+ D(T/q)x c (k + l/q). 0 7 ) 

Substituting (15) and (16) into (17) and rearranging, 
x p (k + 2/q) = A 2 p (T/q)x p (k) 

+ [A p (T/q)D(T/q) + D(T/q)A c (T/q)]x c (k), 
where, in the interest of clarity, it is noted that 

A(T) = A 9 (T/q)=[e AT,q ]‘ l , (19) 


for both the plant and control systems. Continuing a similar 
algebraic process, 

x p (,k + 3/ q) = A p (T / q)x p (k) + 

[Al(Tlq)D(Tlq) + A p (Tlq)D(Tlq)A c (Tlq) ( 20 ) 

+ D(Tlq)Al(Tlq)]x c (k), 
and 

x p (k + 4 / q) = A\ (T / q)x p (k) + 
\Al(T/ q )D(T/ q ) + A 2 p (T/q)D(T/q)A c (T/q) 

+ A p (T/q)D(T/q)A 2 (T/q)+ (2 i) 

D(Tlq)Al(Tlq)]x c (k). 

Then, by induction 


x p (k + l) = Al(Tlq)x p (k) + 

[^A;-\Tlq)]D(Tlq)Ar{Tlq)\x c (k). 


( 22 ) 


B. Method 2: Mapping-Matrix Construction by Continuous 
System Simulation 

The role of D in a discrete-continuous control system 
is to map the contribution of current states of the control 
system across the control decision interval (sample period) 
to future states of the plant. Thus, responses of a 
continuous model of the plant, driven by output from the 
(continuous) function generator portion of the D/C control 
system, can be repetitively propagated across the sample 
period, by dynamic simulation, to generate data from 
which D can be constructed. A block diagram of the 
simulation construct is shown in Fig. 2. 

For an n -dimensional plant and an m -dimensional 
control system, the dimension of D is nxm, and m 

simulation runs are required to construct D . Because 
absolute scale is of no theoretical concern for a linear 
system, the columns of an mxm identity matrix can be 
used to define an appropriate m -dimensional initial state 
vector for the function generator for each of the m 
simulation runs. Initial values of the plant states are set to 
zero so that the contribution to plant response from the 



Finally, by comparison with (12), 


Fig. 2. Diagram of Continuous Simulation 









function generator alone is evaluated. Each simulation run 
produces final values for n elements of the plant state 
vector, and D can be constructed from the resulting data 
set. The basic mathematical relationship supporting the 
computational process is 


x p (T) = D(T)x c ( 0). (25) 

Then, 

ram ) (?) = 

[x' r (T) ! *J(T) | ... | . (26) 


Then, by inspection, comparisons, and substitutions among 
(12) and (30-34), it becomes clear that 



(35) 


where D is as defined by (2). Then, an obvious method 
for computing D is to compute A using any of several 
known standard methods [4], and extract D from the 
partitioned result. 


IV. Implementation 


Because the initial control system states are the columns of 
an m-dimensional identity matrix, 

the convolution integral is given by 

[x\(T ) ; x\(T ) ; ... j *;(r)] ( „, = 

n<«.,(r)[4(0) i *» ! - ! *;(»)],_, • <27) 

where the x‘ P ( T ) are the «-dimensional plant state vectors 

at time t = T which have been propagated by the i‘ h 
execution of the continuous system dynamic simulation of 

Fig. 2 with initial conditions x‘ p (0) and x‘ c (0) given by 

[*i(o) ; jcj(o) ; — | =o, (28) 

and 


Algorithm functionality for Method 1 is defined by (23) 
and (24). Using this definition as a guide, the process was 
programmed directly using the core capability of Matlab®. 

For Method 2, a combination of Matlab® and 
Simulink® was used to build the simulation described in 
Fig. 2, manage its execution, store results, and perform 
subsequent computations as indicated by (27)-(29). 

Algorithm functionality for Method 3 is defined by 
(32), (34), and (35), and was implemented in Matlab®. 
However, two different mathematical processes were used 
to compute the state transition matrix A defined in (35). 
The first is an SPP series-expansion approach derived from 
truncating a series expansion of A analogous to that 
shown in (19), while the second uses the Matlab® library 
function expm ( ) . 


[*c(0) j *c(0) ! •" ! = • (29) 

C. Method 5: Sample-Period-Partitioned Single Series 
Expansion 

A combined state-variable representation of the plant 
and D/C control system function generator depicted in Fig. 
2 can be represented by 

X = A sys X, (30) 


V. Functional Verification 

Initial verification was accomplished by using all three 
methods to evaluate D for three specific cases whose 
closed-form solutions are available [5] and [6]. With 
reasonably proper selection of algorithm parameters, all 
methods appear to function correctly. Algorithm 
parameters include number of terms retained in series 
expansions, number of sample-period partitions, and 
maximum allowable integration step size in the case of 
simulation. The closed-form test cases are outlined below. 



Then, the exact discrete-time representation is 


X((k + \)T) = A iys X(kT), 

1 _ p A sys T 

A sys ~ e 


(31) 


(32) 


(33) 

(34) 


A. Case I: Harmonic Oscillator Plant with Constant-Only 
D/C Control Basis Functions (ZOH control). 

This system consists of an un-damped second-order 
oscillatory plant and a control system whose function 
generator is only capable of producing actuator commands 
that are constant across the sample period. This case is 
equivalent to the ZOH (Zero Order Hold) associated with 
traditional discrete-time control. The matrices that define 
this example system are given in [5] as 


{A p , B p } 




where 


(36) 



and 


and a control function generator defined by 


K , cj ={ o , i }, 


(reference) expression for D in this case is [5] 
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where CO nc is the frequency associated with the control 
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ref 
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Relative computational accuracy can be measured as an 
error norm defined by 


D-D 


error = 


ref\\ 


A 


(39) 


ref 


where | ( )|| can be any suitable matrix norm, and was 
evaluated using the “largest singular value” (default) option 
associated with the Matlab® function norm( ) . Algorithm 

error results, along with pertinent parameter values, are 
shown in Table I. 

B. Case II: Harmonic Oscillator Plant with Linear-in- Time 
(LiT) Control Basis Functions 

This system, like Case I, consists of an un-damped 
second-order oscillatory plant as described by (35), but the 
control system function generator is capable of producing 
LiT actuator commands (constant plus ramp) across the 
sample period. The function generator in this example is 
defined by [6] 


K , Q } = 


il 


[i 


(40) 


system basis functions. The closed-form expression for 
D for this system is [6] 



^W nc T -&\n{o) nc T)} 
-©^[cosK^-l] 

®- 2 r 2 /2+©,> 0 s(© M r)-i1 
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Computational accuracy is measured as before, and 
algorithm parameters and error results are shown in Table 
III. 


TABLE I 

Functional Test Results, Harmonic Oscillator, 
ZOH Control 


Parameter 

Set# 

i 

2 

3 

4 

Plant 

Parameters 

=1 

= 2 

a np =3 

=4 

Control System 
Parameters 

II 

p 

o 

1 — > 

II 

O 

r = i 

II 

O 

Algorithm 

Parameters* 

N = M = \Q 
4=1 

N = M = 4 
5=10 

=M = 10 
5=100 

N = M = \0 
5 = 1000 

Error Norm 

Double Series 

1 .4427e-015 

2.6623e-016 

2.3457e-01 5 

3.6833e-01 5 

Simulation 

1.4878e-015 

7.8487e-017 

2.3969e-01 3 

2.3633e-01 2 

Single Series, 
SPP 

1 .4413e-015 

2.6512e-015 

1 .6612e-01 5 

2.3159e-014 

Single Series, 
exp/77 () 

1 .4544e-01 5 

1 ,4062e-01 6 

7.7714e-01 6 

3.201 9e-01 5 


*M applies only to dbl. series method; max. stepsize=5e-04 for all simulations. 


and the closed-form expression for D from [6] is 


D = 


^[l-cosK/)] 

®«J si 


^K P ^-sinK p 7)]l 
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.(41) 


Computational accuracy is measured according to (39); 
algorithm parameters and error results are shown in Table 


II. 


C. Case III: Double Integrator Plant with Three-State 
Control function Generator. 

This system consists of a double integrator plant, 
defined by 



TABLE II 

Functional Test Results, Harmonic Oscillator, 
LiT Control 


Parameter 

Set# 

i 

2 

3 

4 

Plant 

Parameters 

II 

3^ 


=3 

a np ~ 4 

Control System 
Parameters 

II 

o 

b 

II 

p 

r = i 

o 

II 

Algorithm 

Parameters* 

N=M = 6 

5=1 

N ~M ~1 
5= 10 

N=M = 10 
5= 100 

N=M = 10 
5 = 1000 

Error Norm 

Double Series 

1.4631e-015 

1.3963e-016 

1 .7847e-015 

4.8308e-015 

Simulation 

1.6502e-015 

8.3737e-017 

1.8350e-013 

7.9146e-013 

Single Series, 
SPP 

1 ,4899e-015 

1.3936e-016 

1 .2828e-015 

7.7684e-015 

Single Series, 
expm () 

1.5485e-015 

1.4172e-016 

5.9622e-016 

1 ,4203e-01 5 


\*M applies only to dbl. series method; max. stepsize=5e-04 for all simulations. 





TABLE III 

Functional Test Results, Double Intrgrator, 
3-State Control 


Parameter 

Set# 

i 

2 

3 

4 

Plant 

Parameters 


co =2 

nc 

= 3 

xr 

II 

* 

3 

Control System 
Parameters 

II 

p 

b 

II 

p 

T = 1 

II 

o 


Algorithm 

II 

S! 

II 

< 

II 

S 

II 

JV =M = 10 

N =M =10 

Parameters* 

q=l 

II 

O 

o 

o 

II 

CJ1 

q = 1500 

E 

o 

2 

Double Series 

1.4425e-015 

1.3953e-016 

5.9627e-016 

1.8218e-015 

Simulation 

1.6051e-015 

9.8446e-017 

3.7764e-014 

2.0503e-013 

o 

LU 

Single Series, 
SPP 

1.4617e-015 

1.3939e-016 

6.5836e-016 

8.6898e-015 

Single Series, 
expm( ) 

1.4530e-015 

1.8313e-017 

2.9141e-017 

3.6425e-017 


*M applies only to dbl. series method; max. stepsize=5e-G4 for all simulations. 


® eX P 

uses Matlab library ftmction () to compute the 

matrix exponential seems to produce the most accurate and 
efficient process most of the time. As might be expected, 
the simulation method is generally much less efficient than 
the other methods, and its accuracy noticeably diminishes 
as sample period increases. The double series expansion 
method, while less efficient than the single series method, 
appears to have significant potential for improvement at the 
code level. It also provides a symbolic formula for direct 
evaluation of the convolution in terms of plant, control 
system, and algorithm tuning parameters. It is possible that 
this feature could be of farther use in D/C control 
development. For the constant-coefficient and constant- 
sample-period case, the computation is needed only once, 
and can be done before-the-fact, so efficiency may be of no 
real concern. 

VI. Summary and Conclusions 

Three fundamentally different methods for numerical 
evaluation of the dual-kemal counter-flow convolution 
integral associated with Discrete/ Continuous control theory 
have been derived and implemented. The first method relies 
on series expansion and truncation for the individual plant 
and control system state transition matrices, and also makes 
use of the state transition matrices to partition the sample 
period so that the length of time for which the series 
expansions must be valid can be much less than the sample 
period. The second method requires simulation of the 
system to generate plant time responses across the sample 
period as a function of control system initial conditions. The 
simulation is run as many times as there are control system 
states. A third method relies on a single series expansion 
and truncation for the total system state transition matrix 
followed by extraction of the convolution as a partitioned 
block of the matrix result. This method also incorporates 


sample period partitioning. All three methods appear to 
accomplish the objective in a functional sense; however, 
they are significantly different, and are likely to have 
different strengths and weaknesses with respect to 
accuracy, efficiency and utility in future developments 
associated with Discrete/ Continuous control. 
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